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ps) 'Abstract 

this paper we study the macroscopic conduction properties of large but finite binary networks with conducting bonds. 
fS| 'By taking a combination of a spectral and an averaging based approach we derive asymptotic formulae for the conduction 
5^ 'in terms of the component proportions p and the total number of components N . These formulas correctly identify both 
Onthe percolation limits and also the emergent power law behaviour between the percolation limits and show the interplay 
^between the size of the network and the deviation of the proportion from the critical value of p = 1/2. The results 
compare excellently with a large number of numerical simulations. 
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1. Introduction and summary 

Large but finite binary networks comprising disordered 
"mixtures of two interacting components can arise both di- 
rectly, in electrical circuits P, 0, y, 0] or mechanical struc- 
tures 0, and as models of other systems such as disor- 
idered materials with varying electrical 0, thermal, me- 
"chanical or even geophysical properties in the micro-scale, 
coupled at a meso-scale Q- They are prototypes of many 
forms of complex systems which are often observed to have 
.macroscopic emergent properties which can have emergent 
'power-law behaviour over a wide range of parameter values 
which is different from any power law behaviour of the in- 
dividual elements of the network, and is a consequence of 
the way in which the responses of the components combine. 
"For certain ranges of parameters we see the extensively 
studied percolation type of behaviour Q , in which the over- 
all conductance is directly proportional to the individual 
component conductances with a constant of proportional- 
"ity dependent both on the component proportion and on 
the network size. As well as being important in their own 
right, such large binary networks also provide a useful test 
|bed for identifying different types of emergent behaviour, 
determining what causes it, finding the range of parame- 
ters over which it applies and addressing the fundamental 
question of which aspects of a complex system, such as 
the number and proportion of the components, lead to, 
and influence, the emergent behaviour. In this paper, we 
will combine a spectral analysis, motivated by [9!]), of the 
(partly random) linear operators (Kirchhoff-typc matrices) 
associated with the network, with the averaging methods 
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described in [J], to derive an asymptotic formula for the 
emergent network admittance, that includes both the ef- 
fects of the component proportion p and the network size 
TV. 

As an example, we consider (a set of random realisa- 
tions of) a binary square network comprising a random 
mixture of N conducting bonds which are either chosen 
to have a constant conductance yi = 1/R or a variable 
complex admittance y2 = iuiC , which is directly propor- 
tional to an angular frequency u. If p is the occupation 
probability for choosing a y2 component (and (1 — p) the 
occupation probability for yi), in the limit of large TV or 
for averages over large numbers of systems, it directly de- 
termines the proportion of y2 to be approximately p (and 
yi to be (1 — p)). This network, when subjected to an 
applied alternating voltage of angular frequency uj, has 
macroscopic features, such as the total admittance Y{uj). 
Over a wide range of frequencies < wi < a; < a;2, the 
admittance displays power law emergent characteristics so 
that the magnitude of the complex admittance 1^ | is pro- 
portional to for some a{p). The effects of network size, 
and component proportion, arc important in that uji and 
UJ2 depend upon both p and TV and it is well known Q that 
the casep = 1/2 is a critical value {pc) for two-dimensional 
square networks, lip ^ 1/2 and TV is sufficiently large then 
this problem can be studied by averaging Q , with wi — > 
and a;2 — ^ oo as p — ^ 1/2. In contrast, if p = 1/2, uji is in- 
versely proportional to TV and uj2 directly proportional to 
TV, as TV increases to infinity. For < w < wi and w > a;2 
percolation type behaviour is observed for which Y is pro- 
portional to either yi or y2 with a constant of proportion- 
ality subtly dependent on \p — 1/2 1. Hence we see in this 
system (i) an emergent region with a power law response 
depending on the proportion but not the arrangement or 
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number of the components (ii) a more random region (iii) 
a transition between these two regions at frequency values 
which depends on the number and proportion of compo- 
nents in the system. An illustration of the different types 
of observed response is shown in Figure [T] 

The purpose of this paper is to give insight into this be- 
haviour by obtaining asymptotic formulae for the expected 
response curves. We compare and extend results obtained 
by two complementary methods, one based upon averag- 
ing Q and the other based on properties of the spectrum 
of certain operators i9| . The averaging method works well 
when p 7^ 1/2 and N — > oo, and the spectral method, in 
contrast, works well for the case of p = 1/2 and large, finite 
N. By combining the averaging method with the spectral 
results we will also give approximate asymptotic formulas 
(j68l69|) for Y, valid for general values oi p « 1/2 and 
sufficiently large N, which demonstrates the interplay be- 
tween network size and component proportion. The spec- 
tral method is based both on rigorous results concerning 
the poles and zero distribution of the function Y{uj) and 
on certain semi-empirical results on the regularity of their 
statistical distribution. From these observations we de- 
rive, in the case of p — Pc — 1/2 an asymptotic form for 
the response which gives results almost indistinguishable 
from the numerical simulations, showing that a power law 
response in which is equal to y/uC/R, is observed 

over a range 

1/NCR ^ uji < Lu < L02 = N/CR. 

For ljjCR > N or ljCR < 1/N this is replaced by a 'perco- 
lation' type response for which ji'l is proportional to one 
of 



1/{VNR) OT VNujC for w < 1, 
and 

Vn/R or ujC/Vn for w > 1. (1) 

When p 1/2 we find, both numerically and asymp- 
totically, that for p 1/2 very similar results are obtained 
to those obtained for p = 1/2 if 1 < TV < \l/2-p\^^. For 
larger values of N the response becomes independent of 
N and asymptotic to a p dependent response. For values 
of p not too close to 1/2 this behaviour is approximately 
predicted by the Effective- Medium- Approximation (EMA) 
homogenisation approach based on averaging [l3| , [j| • For 
large and small values of u, the EMA predicts percolation 
limits in which is proportional to 

{l~2p)/R for w<wi, 
and 

l/{{l-2p)R) for UJ-»UJ2 

if p<l/2, (2) 
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Figure 1: (color online) Examples of the type of response observed 
for the proportion p of capacitors (a) below and (b) above the per- 
colation threshold p = pc = 1/2. 



and 



ujC/{2p~l) for a;<a;i, 
and 

{2p — l)wC for uj ^ UJ2 
il p> 1/2. 



(3) 



The percolation limits uji{p) < i02{p) satisfy LOi{p) — > 
0,a;2(p) — > oo as p — > 1/2. In the power law emergent 
region wi < a; < we see power law behaviour propor- 
tional to oj". For C-R networks we will show that the 
EMA implies that 



1 

a = 

2 



2V2v/l -eV2 



' € = l- 2p 



and for R-R networks with y2 = fJ-yi with /i real, we see 
power law behaviour proportional to /i^. As described 
above, we will consider a combination of the EMA and 
spectral approaches which allows for finite size effects and 
gives formulas (j68l69l) for Y involving N and p, in which 
both (jll2l3p arise as special cases. We note, however, that 
the EMA prediction is not particularly good in the limit 
oi p — 1/2. In particular, it has been observed empirically 
, that rather than having percolation limits proportional 
to |l/2 - p\ or |l/2 - in the limit of |l/2 - p| < 1 

they are more closely approximated by expressions of the 
form 

|y|^|l/2-p|±^ /3«1.3. (4) 
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The layout of the remainder of this paper is as fol- 
lows. In Section[2]we will give a series of numerical results 
for a general binary network with admittances yi and 2/2, 
which illustrate the various points made above on the na- 
ture of the network response and will look at both power 
law emergent behaviour and at percolation responses. In 
Section |3] we will formulate the matrix equations describ- 
ing the network and will show how the poles and zeros of 
the admittance function interlace. In Section|3]we will dis- 
cuss, and derive, a series of statistical results concerning 
the distribution of the poles ad zeros. In Section[5]we will 
use these statistical results to derive a precise asymptotic 
form of the admittance y of a general binary network, 
when p = 1/2 and N is large. In Section [5] we review the 
(classical) averaging method for N = 00 which gives an ex- 
cellent estimate when p is not too close to 1/2, and will also 
consider a combination of this method with the spectral 
method for finite N, leading to the formulas (|68l69p for the 
response for all p and sufficiently large N. In Section [7] we 
compare the predictions of the asymptotic formulae with 
numerical computations of the network responses. Finally 
in Section[8]we will draw some conclusions from this work. 

2. Network models and their responses 

In this section we consider basic models for compos- 
ite materials and associated random binary electrical net- 
works with bonds having admittance yi and y2, and present 
the graphs of their responses. In particular we will look in 
detail at the existence of a power law emergent region, and 
will obtain empirical evidence for the effects of network size 
A'^ and capacitor proportion p, on both this region and the 
'percolation behaviour' when = I2/2/1/1I is either large 
or small. 

2.1. Composite materials and their properties 

An initial motivation for studying binary networks comes 
from models of composite materials. Disordered two-phase 
composite materials are found to exhibit power-law scaling 
in their bulk responses over several orders of magnitude in 
the contrast ratio of the components [lllj 01 , and this effect 
has been observed [l|, [l^ in both physical and numerical 
experiment son conductor-dielectric composite materials. 
In the electrical experiments this was previously referred 
to as "Universal Dielectric Response" (UDR), and it has 
been observed fl^ that this is an emergent property 
arising out of the random nature of the mixture. The 
same response is also found in numerical studies using 2D 
lattice structures with bonds randomly assigned to have a 
conductivity which is either constant or linearly variable in 
a parameter. These studies reveal that the emergent scal- 
ing is a property of a large number of complex system that 
can be represented as such a binary percolation network. 
A simple model of such conductor-dielectric mixtures with 
fine structure is a large electrical circuit representation, 
replacing the constituent conducting and dielectric parts 



with a linear C-R network with iV ^ 1 resistors and ca- 
pacitors, respectively forming the bonds in this network. 
Similarly, large mixtures of materials with different resis- 
tance can be modelled by R-R resistor-resistor networks. 
For a binary disordered mixture, the different components 
can then be assigned randomly to the bonds on the lat- 
tice [1^. In most previous studies a 2D square lattice has 
been used, with bonds assigned randomly as either C or 
R, with probability p, 1 — p respectively. The components 
are distributed in a two-dimensional lattice between two 
bus-bars. On of which is grounded and the other is raised 
to a potential V{t) = Vexp{iujt). This leads to a cur- 
rent I(t) = I{oj) exp(za;i) between the bus-bars, and we 
measure the macroscopic (complex) admittance given by 

Y{uj) ^ I{lj)/V. 

This approach is closely related to percolation models and 
a large review of this and binary disordered networks can 
be found in d, 0, |1] . There are many advantages to using 
C-R network representations of these types of system. In 
particular, widely available circuit simulation software can 
be used, which makes use of the available efhcient sparse- 
matrix techniques in solving the equations of the system. 
Additionally, algorithms based on the Frank-Lobb reduc- 
tion techniques can be used on the 2D square lattice 
representations to solve large systems efficiently [1} . These 
techniques were used in various studies to show that the 
PLER exists in an y b inary random network with variable 
contrast ratio 13, S S [3 • This allows many differ- 
ent simulations to be made of different realisations of the 
circuit with randomly assigned resistors and capacitors. 
Furthermore, finite element calculations reported in [l^ 
indicate that the response of the full material is very close 
to that of the network model of that material. 

2.2. Percolation and power-law emergent behaviour 

When a; ^ 1, the capacitors act as open circuits and 
conduction occurs predominantly through the resistors, 
having far higher admittance than the capacitors. The 
circuit then becomes a percolation network in which the 
bonds are either conducting with probability (1— p) or non- 
conducting with probability p. The network only conducts 
only if there is a percolation path from one electrode to 
the other, and it is well known [19] that, for 2D square lat- 
tices the critical percolation probability Pc — 1/2. Thus, 
in cases where the probability of the non-conducting phase 
p > Pc — 1/2, the conducting phase has a very low proba- 
bility for such a percolation path to exist. 

In contrast, if p < 1/2 then such a path exists with 
probability approaching one as the network size increases. 
The case of p = 1/2 is critical with a 50% probability that 
such a path exists. This implies that if p < 1/2 there is 
almost certainly only a resistive conduction path and for 
low frequencies the overall admittance is independent of 
uj. In contrast, if p > 1/2 then there is almost certainly 
no path through the resistors and the admittance is ca- 
pacitance dominated and directly proportional to cj. If 
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p = 1/2 then half of the reahsations will give an admit- 
tance response independent of lo and half an admittance 
response proportional to lo. When w 3> 1, we see an op- 
posite response. In this case the capacitors act as almost 
short circuits with far higher admittance than the resis- 
tors. Again we see percolation behaviour with the resis- 
tors behaving approximately as open circuits in this case. 
Thus \i p > 1/2 any conduction path is most likely solely 
capacitative, with the resulting overall admittance being 
proportional to w, and if p < 1/2 a response independent 
of ijj. The case of p = 1/2 again leads to both types of 
response having equal likelihood of occurrence, depending 
upon the network configuration. Note that this implies 
that if p = 1/2 then there are Jour possible qualitatively 
different types of response for any random realisation of 
the system. For intermediate values of w the values of the 
admittance of the resistors and the capacitors are much 
closer to each other. It is here that we see power-law 
emergent behaviour (PLER). This is characterised by two 
features: (i) an admittance response |y| that is propor- 
tional to a;",a « p over a range uj G {uji,ijJ2) and which 
displays a strong symmetry in the behaviour for small and 
large values oi uj. (ii) when p — 1/2 a response that is 
not randomly dependent upon the network configuration. 
Figures [2] (a) and (b) plot the admittance response as a 
function of lj in the cases oi p — 0.4, p = 0.6 and Figure [3] 
for p = 1/2. The figures clearly demonstrate the forms of 
behaviour described above. Observe that in all cases we 
see quite a sharp transition between the percolation type 
behaviour and the PLER behaviour as uj varies. 

We have seen above how the response of the network 
depends strongly upon p. It also depends upon the net- 
work size and this effect is critical if p = 1/2. Fig- 
ure 0] shows the response for the critical value p — 1/2 
for different values of N . Observe that in this case the 
width of the power-law emergent region increases appar- 
ently without bound, as N increases, as do the magnitude 
of the responses for small and large frequencies. From 
these graphs, it is apparent that in this critical case the 
upper limit or the PLER is proportional to N and the 
lower limit proportional to 1/iV. We can very roughly mo- 
tivate the result for p = 1/2 as follows. Suppose that oj 
is small so that the capacitors essentially act as open cir- 
cuits. Imagine for a single percolation path through all 
of these capacitors comprising a chain of resistors, then 
this will have an approximate length of \/N resistors and 
hence a conductance of 1/{\/NR). In contrast, if there is 
a dual path of capacitors going from top to bottom of the 
network, interrupting the resistors, then each resitive path 
has conductance iwC and there are ^/N of these in paral- 
lel, so that the overall conductance is ^/NiwC. In Figure 
[5] we plot the response for p = 0.4 and again increase N. 
In contrast to the former case, away from p — 1/2, the size 
of the power-law emergent region appears to scale with N 
for small N before becoming asymptotic to a finite value 
for larger values of N , consistent with formulae (|2l3p . 
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Figure 2: (color online) Typical responses of network simulations for 
values of p 7^ 1/2 which give qualitatively different behaviour so that 
in the percolation region with oj <^ 1 or (.j 3> 1, we see resistive 
behaviour in case (a) and capacitative behaviour in case (b). The 
figures presented are density plots of 100 random realisations for a 
20 X 20 network. Note that all of the realisations give very similar 
results. 
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Figure 3: (color online) Responses for 100 realisations at p = 1/2 
showing four different qualitative types of response for different re- 
alisations. Here, about half of the responses have a resistive perco- 
lation path and half have a capacitive one at low frequencies with 
a similar behaviour at high frequencies. The responses at high and 
low oj indicate which of these cases exist for a particular realisation. 
The power-law emergent region can also be seen in which the admit- 
tance scales as y/ui and all of the responses of the different network 
realisations coincide 
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Figure 4: (color online) The effect of network size on the width 
of the power-law emergent region for p = 1/2, in which we see this 
region increasing without bound as N increases. 
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Figure 5: (color online) The effect of the network size TV on the 
power-law emergent region for p = 0.4, in which we see this region 
becoming asymptotic to a finite set a,s N —> oo. 
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2.3. The effects of network size and capacitor proportion 
To compare these results and to investigate the inter- 
play between network size and the proportion p, we con- 
sider for p < 1/2, the dynamic range of the response for 
those realisations which have a resistive percolation path 
for both low and high frequencies (that is with probability 
one if p < 1/2 and probability 1/4 if p = 1/2). We define 
the dynamic range Y{N,p) by 



Y 



\Yl 



Yl 



|Y|(cj ^ oo) 



In Figure|n](a) we plot F as a function of N for a variety of 
values oi p < 1/2. We see from this figure that if p = 1/2 
then Y is directly proportional to TV for all values of N. 
In contrast, Hp < 1/2 then Y is directly proportional to 
N for smaller values of N and then becomes asymptotic 
to a finite value Y{p) as N — > oo, with the asymptotic 
behaviour occurring when N > (1/2 — The for- 

mulae ([BSinnilj derived in the final section by combining the 
EMA and spectral results, imply that Y is approximated 
by where /3 satisfies the quadratic equation 



N 



(1 - 2p)/3 -1=0. 



(5) 



This gives reasonable qualitative agreement with the 
calculations presented in Figure [5] with Y ^ N ioi smaller 
values of N and Y ^ Y{p) as iV — >■ cxj. However, we do 
have to exercise a certain degree of caution in applying 
this formula. In Figure |6] (b) we present Y{N,p) as a 
function of p as p ^ 1/2, showing the limiting value Y{p) 
of Y{N,p) as N is increased to infinity. We see in this 
figure that whilst the estimate Y{p) ^ {1 — 2p)~^ is fairly 
accurate, a much better estimate in the limit of p — > 1/2 
is given by 

Yip) ^ (1 - 2p)-2-6 

which is consistent with known empirical results on perco- 
lation 

3. Linear circuit analysis 

We now describe in detail how the disordered material 
is modelled by a general electrical network model with 
two types of bond of admittance yi and 2/2 in respective 
proportions 1 — p and p. These will have admittance ratio 
/X = 2/2/2/1- For a resistor-resistor (R-R) network with yi — 
1/Ri and 2/2 = I/-R2 we have 



fi = R1/R2 is real and positive. 



(6) 



For a capacitor-resistor (C-R) network with yi = 1/R and 
2/2 = iujC we have 



/i = iujCR is purely imaginary. 



(7) 



For a capacitor-inductor (C-L) network /i is a real negative 
number but we do not consider this case here. Our interest 
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Figure 6: (color online) (a) Variation of the dynamic range Y = 
\Y\max /\Y\min as a function of N and p, (b) The value of Y(p) as 
a function oi p ^ 1/2 comparing the estimates (1 — 2p)~^ and (1 — 
2p^— 2.6 Pqj. gach value of p the vertical sequence of dots represents 
calculations of Y for increasing values of N . 



will be in how the overall admittance of the system varies 
as the admittance ratio itself varies, and considering how 
this can be determined in terms of the poles and zeros of 
the admittance function Y{^). 

3.1. Linear circuit formulation 

Now consider the 2D N node square lattice network 
shown in Figure [71 with all of the nodes on the left-hand- 
side connected via a bus-bar to a time varying voltage 
V{t) = T/e*"* and on the right-hand-side via a bus-bar 
to earth iV). We assign a voltage Vi with i — 1 ... to 
each (interior) node, and set v = {vi,V2tV^ . . .vn^^ ■ We 
also assume that adjacent nodes are connected by a bond 
of admittance z/ij- S {2/1,2/2}- The current from the node 
i to an adjoining node at j is then given by lij where 
li.j — {vi — Vj)yij. From Kirchhoff's current law, at any 
interior node all currents must sum to zero, so that 



0. 



(8) 



If i is a node adjacent to the left boundary then certain 
of the terms Vj in will take the value of the (known) 
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Figure 7: (color online) Illustration of an example Resistor Capacitor 
circuit. 



applied voltage V{t). Similarly, if a node is adjacent to the 
right hand boundary then certain of the terms Vj in 
will take the value of the ground voltage 0. Combining all 
of these equations together leads to a system of the form 



= V{t)h = Fe*"*b, 



(9) 



where K = K{lj) is the (constant in time) N x N sparse 
symmetric Kirchhoff matrix for the system and the adja- 
cency vector b = b(a;) is the vector of the admittances 
of the bonds between the left hand boundary and those 
nodes which connected to this boundary, with zero entries 
for all other nodes. As this is a linear system^ we can take 
V = Ve*'^* so that the (constant in time) vector V satisfies 
the linear algebraic equation 



KV = Vh. 



(10) 



If we consider the total current flow / from the LHS 
boundary to the RHS boundary then we have 

/ = h^{Ve -V) = Vc- b'^V, 

where e is the vector comprising ones for those nodes 
adjacent to the left boundary and zeroes otherwise and 
c = b-^e. Combining these expressions, the equations de- 
scribing the system are then given by 



KV^hV = 0, cV-h^V = I. 



(11) 



The bulk admittance Y{fi) of the whole system is then given 
hy Y = I /V so that 



Y{fi) ^c-h^K-'h 



(12) 



Significantly, the symmetric Kirchhoff matrix K can be 
separated into the two sparse symmetric N x N component 
matrices K — Ki + K2 which correspond to the conduc- 
tance paths along the bonds occupied by each of the two 
types of components. Furthermore, using /i — 2/2/2/1, we 
have 

Ki ^ 2/iLi and K2 = 2/2^2 = A*2/i-^2 (13) 



and hence 

K = yiLi + ^2/1^2, 

where the terms of the sparse symmetric connectivity ma- 
trices Li and L2 are constant and take the values 1, 0, —1. 
Note that K is a linear ajfine function of /i. Furthermore, 

A = L1+L2 

is the discrete, positive definite symmetric, negative Lapla- 
cian for a 2D lattice. Similarly we can decompose the ad- 
jacency vector into two components bi and b2 so that 



bi 



2/iei + 2/262 = 2/iei + ^2/1^2, 



where ei and 62 are orthogonal vectors comprising ones 
and zeros only corresponding to the two bond types ad- 
jacent to the LHS boundary. Observe again that b is a 
linear affine function of fi. A similar decomposition can be 
applied to the scalar c = yici + fJ.yiC2. 

3.2. Poles and zeros 

To derive formula for the expected admittances in terms 
of the admittance ratio /i we now examine the structure 
of the admittance function Y{^). As the matrix K, the 
adjacency vector b and the scalar c are all affine functions 
of the parameter fjL it follows immediately from (fT2|) and 
Cramer's rule applied to (ITT]) that the admittance of the 
network K(/i) is rational function of the parameter /x, tak- 
ing the form of the ratio of two complex polynomials P(^) 
and Q{ii) of respective degrees r < N and s < A^, so that 



Q{fi) qo + qi^ + q2^^^ + ... qrfJ.'^ 



(14) 



Pip) po +PlM +P2M^ + • • -PsM* 

We require that po 7^ so that the response is physically 
realisable, with Y{^) bounded as w and hence ^ — )■ 0. Sev- 
eral properties of the network can be immediately deduced 
from this formula. For convenience, we look at a C-R 
network, although similar results arise for R-R networks. 
First consider the case of uj small, so that fi = iuCR is 
also small. From the discussions in Section [2j we predict 
that either there is (a) a resistive percolation path in which 
case Y{^) ~ /i*^ as ^ — >■ or (b) such a path does not exist, 
so that the conduction is capacitative with Y{fj.) ^ /i as 
/i — > 0. The case (a) arises when 7^ and the case (b) 
when qo = 0. Observe that this implies that the absence of 
a resistive percolation path as — > is equivalent to the 
polynomial Qip) having a zero when fi = 0. Next consider 
the case of a; and hence /i large. In this case 
qr 



Ps 



as /i — ^ 00. 



This time we may have (c) no capacitive path at high 
frequency with response Y{fi) ~ /x" as /i — 00, or the 
existence of a capacitative path with Y{fi) ^ fi. In case 
(c) we have s = r and Pr ^ and in case (d) we have 
s = r — 1 so that we can think of taking pr = 0. Accord- 
ingly, we identify four types of network defined in terms of 
the percolation paths for low and high frequencies, which 
correspond to the cases (a), (b), (c), (d) so that 
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(a) 




(b) 


Po = 


(c) 




(d) 


Pr ^0 



Both the polynomials and can be factorised 

by determining their respective roots fJ.p^k, k — 1 . . . s and 
l^z,k, k = 1 . . .r which are the poles and zeroes of Y(fi). 
We will collectively call these poles and zeroes the reso- 
nances of the network. Our analysis of the network will 
rely on determining certain statistical and other proper- 
ties of these resonances. Note that in Case (b) we have 
^z,i = 0. Accordingly the network admittance can be ex- 
pressed as 



Y{^i,N)^DiN) 



n (a* - 

k=l 

s 

n (/^ - 

fe=i 



(15) 



Here D{N) is a function which does not depend on fi but 
does depend on the characteristics of the network. 

3.3. Location of the resonances 

We now proceed to prove some rigorous results con- 
cerning the location of the poles and zeroes. We firstly 
note that the number r of poles/zeroes can be substan- 
tially less than N. This is due to the formation of clus- 
ters of components in the lattice which are isolated from 
the boundaries 0. Such component clusters lead to reso- 
nances at infinity or zero, depending on which component 
the clusters are made of. Comparing (|12p and (ITSI) . it can 
be seen immediately that the poles are precisely the roots 
of the determinant of the Kirchhoff matrix K. This matrix 
has the form K = Ki -\- K2 = yi{Li -f with Li, L2 
constant and symmetric (though not necessarily positive- 
definite) and Li + L2 = A. The poles are then given by 
— 1 times the eigenvalues of the matrix pencil (Li, L2), so 
that the values Hp^k, and the corresponding vectors Vp^k, 
satisfy the linear equation 

(Li -I- fJ.p^kL2)^p,k = with Wp^k ^ 0. (16) 

As Li + L2 = A this then implies that 

(Li(l - [ip^k) + fJ.p,k^)vp,k = 



so that 



(Li + iip^k^/ (1 - tipM))^p.k = 0. 



It follows immediately from the symmetry of Li and 
the fact that A is a symmetric positive definite operator, 
that /Xp_fe/(1 — fip^k) is real. The negativity of /ip^fc follows 
from the fact that the network has a bounded response. 

In contrast, using (|12l) and PT|) . the zeros are those 
values of /i — ^z,k, with corresponding vectors Vz.k which 
satisfy the simultaneous equations 



(il + Aiz,fci2)Vz,A; = b, b'^V, 



z.k 



c, 



so that the current / is zero. The condition for the zero 
can thus be put into an extended matrix equation of the 
form 

r , 1 

0. 

As if, b and c are all affine functions of /i, this leads to the 
corresponding eigenvalue problem for the zeros /i^ k given 
by 



" K ~h 






-b^ c 




1 



yi{Li -\- fiz^kL2) 2/i(-ei - ^^,fce2) 
yi(-e|' - Aiz,fce^) yiici + Hz.kC2) 



v.,fc=0, (17) 



with Vz^k 7^ 0, so that the zeros are the eigenvalues ^z,k 
of this extended matrix pencil. Both problems (|16l) and 
PT)) can be transformed into standard eigenvalue prob- 
lems. From the previous reasoning, problem (1161) is equiv- 
alent to the problem 

l^p,k 



p,k 



l^p,k 



1 



-Av, 



p,k 



(18) 



and using the fact that A is a symmetric positive definite 
matrix, the Cholesky decomposition [20] A = LL^ exists, 
where L is a lower triangular matrix. Thus we can rewrite 
(dH]) as a standard eigenvalue problem 

L^^LiL^'^Wp^k ^ Cp.k^P,k, with Wp^fe 7^ 0, (19) 



Note that as 



where Wp,fc = L^^p^k and Cp,fe = ^^'^ . 

lJ'P,k — ^ 

fip^k < it follows that < (p^k < 1- Similarly, (flT)) can 
be written as 



-ei 

Cl 



IJ'z.k 



l^z.k 



1 



-ei - 62 

Cl + C2 



Using the Cholesky decomposition of A again we can de- 
fine an extended generalised Cholesky decomposition of 
the extended matrix system by 



Li +L2 



-e^ — e 



(ei +e2)'L 



-ei - 62 

2 Cl + C2 

L 0' 







a 



- 62 



where = ci -I-C2 — (ei + e2)^(LL-^)^^(ei +62) and L is 
a lower triangular matrix. Using this extended Cholesky 
decomposition we rewrite ([TT]) as 



Wz,fe = Cz.k^z.k, with ^z.k ^ 0, 

(20) 

= {L-^LiL-^{ei + 



l^z,k 



lJ'z,k - 1 

e2)-L~-^ei))/a and c = {{61+62)'^ L^^LiL^'^ {ei+ 62) - 
ejL~'^{ei -f ea) - (ei + e2)^L-^ei + ci))/a^. 

Now, by the Cauchy interlacing theorem [21I Theo- 
rem 10.1.1] (see also [1^), the eigenvalues of (fT9)) interlace 
those of ([^, that is 



< Cz,l < Cp,l < Cz,2 < Cp,2 < . . . < Cp,s < Cz,is+1) < 1- 



Furthermore the eigenvalues Cz,k,i = 1, . . . ,s + 1 corre- 
sponding to the zeros are given by the zeros of the function 



A 



Equivalently, the poles and zeros of Y{fi) are all nega- 
tive real numbers, and interlace so that 



(21) 



This result immediately leads to two different interpre- 
tations in the case of an R-R and a C-R network. In the 
case of an R-R network with conductance ratio /i > the 
poles and zeros occur along the negative real axis so that 
we can take fip^k — —Mp_k < etc. Thus, as fi varies 
through positive real values 



nLi(M + M,,fc) 



(22) 



with the values M| > and M|? > 0. For the C-R network. 



jCR, and we can consider y to be a function of w. 



The poles ujp^k of Y{u;) then satisfy iCRujp^k = —Mp^k so 
that they lie along the positive imaginary axis, ditto the 
zeros. Thus, as uj varies through real values then 



r(w) = D{N) 



nLi(- 



(23) 



with the values > and Wj^ > 0. We note that neither 
of the expressions (^51 become unbounded as uj varies 
through real values or as varies through positive real 
values. This is in contrast to the case of an C — L network 
in which the resonances can be real and positive can lead 
to unbounded responses as uj varies. In contrast, we see 
in the C — R and R — R networks, an averaging effect in 
the product terms in these expressions, which leads to the 
emergent behaviours observed in practice. 

4. The distribution of the resonances 

We now look at the distribution of the poles and zeros, 
and draw certain conclusions about their statistical regu- 
larity which will allow us to then compute the asymptotic 
form of the system response. The statistics of the reso- 
nances are most regular in the critical case of p = 1/2, 
allowing us to make very precise estimates of the over- 
all system behaviour in this case, precisely complementing 
the averaging methods which work best when p ^ 1/2. To 
perform these calculations, we note that if we consider the 
elements of the network to be assigned randomly, with the 
components taking each of the two possible values with 
probabilities p and (1 — p), then we can consider the reso- 
nances to be random variables. The poles and associated 



eigenvectors are given by the solutions of the matrix pencil 
equation (see p^ ) 



(24) 



Each realisation of the network, with bonds chosen from 
a Binomial [p, (1 — p)] distribution will give a different set 
of values for ^p^k = —iMp^k = iCRWp^k and we can then 
consider the statistics of this set. We ask the following 
questions: (1) What is the statistical distribution of fip^k 
if N is large? (2) What is the statistical distribution of the 
location of a zero between its two adjacent poles? (3) How 
do fip^i and fip^N- In each case we will find good numer- 
ical evidence for strong statistical regularity of the poles 
(especially in the case oi p = 1/2), leading to answers to 
each of the above questions. 

4-.1. Preliminary observations on the pole locations 

To motivate our answers we start by considering the 
special case oi p = 1/2. In this case the two matrices 
Li and L2 representing the connectivity of the two com- 
ponents have a statistical duality. Statistically, any real- 
isation which leads to a particular matrix ii is equally 
likely to lead to the same matrix L2- Because of this, if 
/X is an observed eigenvalue of the pair (ii, L2) then it is 
equally likely for there to be an observed eigenvalue of the 
pair (L2,Li) with the same eigenvector and with eigen- 
value being precisely 1/fi. Thus in any set of realisations 
of the system we expect to see the eigenvalues ^ and I//1 
occurring with equal likelihood. It follows from this simple 
observation that the variable log(/x) should be expected to 
be a random variable with a symmetric probability dis- 
tribution and with mean zero. It is therefore natural to 
expect that for a large number of realisations, the vari- 
ables \og{Mk,p) should follow a normal distribution with 
mean zero (so that Mp^k has a lognormal distribution cen- 
tred on M = 1). Similarly, if Mp^i is the smallest value 
of Mp^k and Mp^M the largest value then Mp^i = 1/AIp^N. 
In fact we will find that in this case oi p = 1/2 we have 
Mp^i ^ 1/N and Mp^N ~ N. It follows similarly that 
log{Wp^k) is expected to have a mean value of — log(Ci?). 
Following this initial discussion, we now consider some nu- 
merical computations of the distribution of the poles in a 
C-R network for which CR — 10~^. As a first computation 
we consider many random realisations of networks gener- 
ated with a large enough size (typically N ~ 380) to ensure 
good statistics per network. We define S as the number 
of horizontal components in one row of the network; giv- 
ing horizontal and (5* — 1)^ vertical components. The 
number of internal nodes (i.e. not including the boundary 
nodes), which is equal to the dimension of the matrix K, 
is therefore N — S{S — 1); giving the maximum possible 
number of eigenvalues fii. The results of the computations 
are presented in Fig. [8] in which we give a histogram of 
the distribution of the poles Wp^k (on a log-scale in the 
frequency domain) over 100 different realisations of each 
network. These figures clearly indicate that the location 
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of the poles does indeed possess a strong statistical reg- 
ularity, conforming approximately to a log-normal distri- 
bution with mean log(l /Ci?) in all cases. Evidence for 
this is given by comparing the resulting curve with the 
standard Normal distribution with an appropriately cho- 
sen value for the variance. The fitted curves in Figure |8] 
show that the results are close to log-normal for any choice 
of p (provided that N is chosen sufficiently large). When 
the results of the realisations considered above are fitted 
to a log-Normal distribution with probability density func- 
tion P{W) = aexp {-{W - E{W})^/2a^) we find the re- 
markable result that the standard deviation a appears to 
be largely independent of the value of N and to display a 
simple functional relation on p. In Figure ^ we show the 
standard deviation as a function of a for a range of values 
of N and see a good fit to the curve a — ap{l — p) over 
many values of N . As a second calculation we take a single 
realisation of a network with N w 380 nodes and p = 1/2 
and determine the location of the imaginary part of poles 
Wp,k- A plot of the logarithm of the poles, ordered in in- 
creasing size, as a function of k is given in Figure [TUl Two 
features of this figure are immediately obvious. Firstly, 
the terms Wp,k appear to be point values of a regular 
function f{k). Secondly, \og{CRWp_k) shows a strong de- 
gree of symmetry about zero, so that \i 1 < k < N then 
\og{CRWp^k) = if fc = iV/2. Motivated by the discus- 
sion above, we compare the form of this graph with that of 
the error function^ that is we compare erf(log(Ci?T4^p^fc)) 
with 2k/N — 1. The correspondence is very good, strongly 
indicating that log(/) takes the form of the inverse error 
function with an appropriate constant of proportionality. 

4-2. Pole-zero spacing 

As a next calculation we consider the statistical distri- 
bution of the location of the interlacing zeros with respect 
to the poles. In particular we consider the variable r]k, 
which depends on the proportion p given by 



Vk{p) 



logAfp^fc+i - logAf^^fe _ \ogWp^k+i - logVF^^fe 



log Mp.fe+i - log Mp,fc log Wp^k-\-i - log Wp^k 

(25) 

We now establish three symmetry results for the mean 
value fjkip) of rjk, taken over many realisations. 

First symmetry: Assume the zeros are Wz.k, k G [0,N] 
and the poles are Wp^k, k G Define, 



\ogWp^k+i - logtVj 



z.k 



\ogWp^k+i - logWp^k ' 



where. 



So, 



\0gWz,k = -logT4^2,Ar_fc, 
logWp,fe = -logT4^p,Ar_fe+i. 



f]N-k{p) 



^ logW^.fc - log Wp,fc 
~ log VFp,fe+i - log Wp,A 



(26) 



(27) 
(28) 



(29) 




Figure 8: (color online) Distribution of Wp fe values for (a) p = 0.25 
and (b) p = 0.5 over a range of system size: (i) N = 90, (ii) N = 380. 
The curves are fitted to Normal distributions (on a log-scale). The 
variance appears to depend on p but be largely independent of the 
value of N . 
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Figure 9: (color online) The standard deviation tr as a function of p 
showing (T ap{l — p). 



size A' 
1000 



then, 



Hence, 



VN-kip) +rik{p) = 1. 



(30) 



It follows from simple symmetry considerations that Li{p) 
has the same form on average as 2^2(1 — p) and vice versa. 
Hence, if ^,p^k is observed for one realisation with y2 in 
proportion p, then l//ip,fc will be observed when the pro- 
portion of j/2 is 1 — P- A similar result holds for the zeros. 
Noting that. 



f]N-k{v) = VN-k{l-p), 



(31) 




(a) 



-2 -1.5 -1 -0.5 0.5 1 1.5 2 
logio{CRWk)/1.15 



o 



o 




(b) 



Figure 10: (color online) The location of the logarithm of the poles 
as a function of k, for a single realisation of the network, and a 
comparison with the inverse error function. 



Vkip) + r]N~k{i -p) = 1. 



(32) 



Second symmetry: As a second observation we invoke 
duality results due to Keller [l^ (see also Isl ) , in which the 
admittance of a network is compared with that of the dual 
network, in which every bond of the original network is re- 
placed with an orthogonal bond for the dual. Significantly, 
square binary networks are self-dual. A consequence of the 
duality results is that if 



Y{yi,y2) Y{y2,yi) = yi 2/2- 



(33) 



In the case where p = 0.5 it is often argued that as it 
claimed that as F (2/1,1/2) = Y{y2,yi) then Y = ^yi,?/2- 
As we have seen in Section [5J this is only correct in the 
PLER (where the response is unique). It this is not quite 
correct in the percolation region where Y can take one of 
two forms, but it does correctly predict the duality be- 
tween these two forms. It follows from (l33l) that 



D{N) 



nLi(M+M,,fc) ^ 2/12/2 nLi(i/M+Mp,fc) 
nfc=i(A^+^p,fe) D{N) nLi(VM+^^,fc)' 



This can only be true for all ^ if we have the symmetry 
result (taking the ordering of the poles and zeros into ac- 
count) given by 



It immediately follows from ((25|) that asymptotically we 
have the second symmetry 

Vk{p) ='nN-k{p)- (34) 

Third symmetry: Combining (I32p and (1341) . we derive the 
third symmetry 



Vkip) +Vkil-p) = 1- 



In particular, this gives 



r7fe(l/2) = 1/2. 



(35) 



(36) 



The distribution of fjk{p) (over 100 realisations of a C-R 
network) plotted as a function the location of \og{Wp,k) 
for p = 0.3, 0.5, 0.7, is shown in Figure [TT] together with a 
graph of ?7fe(0.3)-|-77fe(0.7). The figures (a),(b) and (c) show 
clearly the refiectional symmetry about the mid-point im- 
plied by dMl). The figure in part (b) (with p = 1/2) is 
particularly remarkable, clearly indicating, as predicted by 
that 77fe(l/2) is equal to 1/2 almost independently of 
the value of log(Wp_fe). There is some deviation from this 
value at the high and low ends of the range due to slower 
convergence to the mean. As well as this there is some 
evidence for a small asymmetry in the results, but the con- 
stancy of the mean near to 1/2 is very convincing. The fig- 
ure in part (d) forp = 0.3 andp — 0.7 clearly illustrates the 
symmetry relation (j35p . We note, however, that for p ^ 
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Figure 11: (color online) Figures showing how the mean value fjf; of 
r]k taken over many realisations of the network, varies with the mean 
value of Wp^k- The four examples show results for (a) p = 0.3, (b) 
p = 0.5, (c) p = 0.7 and (d) 77(0.3) + f?(0.7). 



1/2 the value of fjk of rjk varies with log{Wp,k) in a sym- 
metric distribution (as predicted by ([M)) ) which depends 
approximately quadratically on the value of \og(Wp^k) ■ If 
p > 1/2, fjk takes a value a little less than p in the centre of 
the range when Wp^k — Wmid — l/CR, and a bit greater 
than p towards the ends of the range. The distribution 
is reversed when p < 1/2, as can be seen by comparing 
Figures [TT] (a) and (c), and this is a consequence of 

4.3. Limits of the resonance distributions 

As a final calculation, we consider the number N' of 
the finite non-zero resonances in this case of a C-R net- 
work, and the location of the first non-zero pole and zero 
Wz.i,Wp^i and the last finite pole and zero Wp^N' ,Wz,n' ■ 
As discussed, in the case of p = 1/2 we expect a sym- 
metrical relation so that CRWp^i and CRWp^N' might be 
expected to take reciprocal values. We consider two calcu- 
lations, firstly determining N' /N for a range of values of 
N and of p and secondly calculating the functional depen- 
dence of Wp^i and Wp^N' upon N and p. The value of N' 
can be considered statistically and represents probability 
of a node contributing to the current paths. If we take 
z = N' /N as a function of p for a range of values of N the 
resulting distribution is plotted in Figure [121 We see that 
the shape of this curve is parabolic in p with a maximum 
value for z « 0.8 given when p = 1/2. Indeed, statisti- 
cal arguments presented in Q indicate that the maximum 
value at p = 1/2 is given by A^' = 3 (2 - v^) = 0.804 .... 
We next consider the values of Wp^i and of Wp^N' which 




Figure 12: (color online) The value of N' /N for varying values of p 
and TV. 

will mark the transition between emergent type behaviour 
and percolation type behaviour. A log-log plot of the val- 
ues of Wz,i, Wp^i and of WzM', Wp^N' as functions of A^ 
for the case of p = 1/2 is given in Figure [T51 There is very 
clear evidence from these plots that each of W^p,i and 
Wz,N'i Wp^Ni both have a strong linear dependence upon 
A^ and 1/A^ for all values of N. Indeed we conclude from 
this figure that the following reciprocal relations hold 



CRWzu CR W, 



A^" 
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and 

CR W:,,N', CR Wp^N' - N, 

with an identical scaling for Af^^i, Mp_i, Mz,n' , Mp^^' ■ 
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Figure 13: (color online) Largest W^j^iyWpjf/ and smallest 
Wz,\, Wp^\ zeros and poles for a C-R network with p = 0.5 showing 
linear dependence on A'^ and 1/A'^. 



4- .4^. Summary 

The main conclusions of this section are that there is a 
strong statistical regularity in the location of the poles and 
the zeros of the admittance function. In particular we may 
make the following conclusions based on the calculations 
reported in this section. 



1, 



Mp^k ^ /(k) for an appropriate continuous function 
f{k) where / depends upon p strongly and upon TV 
very weakly. 

fjk{l/2) 1/2 for all values of k. 



Up = 1/2 and if M^,i 



0, then Mp,i,M^ 



5. Asymptotic analysis using spectral results 

We now use the results in the summary of the previous 
section to derive the form of the conductance in the two 
cases of a C-R network and an R-R network. The results 
are sharp for the case of p = 1/2 and N < oo allowing us to 
deduce the asymptotic responses in this case. The formulae 
that we derive will take one of four forms, depending upon 
the nature of the percolation paths. 

5.1. Derivation of the response for general ^ 

We consider the formulae for the value of the admit- 
tance of the binary network 



n 0^ - M^.fc) 
Y{^^)^D{N)^^^ , 

n (a* - Aip.fc) 
fe=i 



(37) 



where the results of section ([3]) imply that ^z,k = — 
, and < M^a < Mp^ < M^,2 < Mp^2 < ... < MpJ{< 
^z(s+i))- Crucially, as /i is either positive or purely imag- 
inary, y is a bounded function for all /i. Here we assume 
that we have s — N' poles, but consider situations with dif- 
ferent percolation responses for |^| large or small, depend- 
ing upon whether the first zero Mz^i = and on the exis- 
tence or not of a final zero (tv'+i)- These four cases lead 
to four functional forms for the conductance, all of which 
are realisable in the case of p = 1/2. In this section we de- 
rive each of these four forms from some simple asymptotic 
arguments. At this stage the constant D{N) is undeter- 
mined, but we will be able to deduce its value from our 
subsequent analysis. Although simple, these arguments 
lead to remarkably accurate formula when p — 1/2, when 
compared with the numerical calculations, that predict not 
only the PLER but also the limits of this region. To obtain 
an asymptotic formula from (j37l) we assume that s = iV' is 
large, and that there is a high density of poles and zeros. 
From the results in the previous section we know that, 
asymptotically, the poles at —Mp^k follow a regular dis- 
tribution and that the the zeroes have a regular spacing 
between the poles, especially in the case of p = 1/2. The 
conclusions of the previous section on the distribution of 
the poles and the zeros leads to the following formulae: 

Mp,k ~ fik), 
Mp^k+i — Afz,fc+i _ ^ 

M,,fc+i ^ f{k) + (1 - (38) 

Here, as we have seen, the function log(/(A;)) is given by 
the inverse of the error function, but its precise form does 
not matter too much for the next calculation. To do this 
we firstly consider the contributions to the product in (|37)l 
which arise from the terms from the first pole to the final 
zero, so that we consider the following product: 



N' 



P^DiN)l[ 



k=l 



(39) 



Note that this product has implicitly assumed the exis- 
tence of a final zero (jv'+i). This is specific to to the 
physical case where there is a percolation path through the 
1/2 bonds but no percolation path through the yi bonds. 
This contribution will be corrected in cases for which such 
a final zero does not exist. Using the results in ([55]) . in 
particular on the mean spacing of the zeros between the 
poles, we may express P as 



N' 

P - D{N) II 

k=l 
N' 



Ai + (/(fc) + (l-4)/'(fc)) 



= m) n 1 



k=l 



^i + f{k) 

(1 - 4)r(fc) 

^i + f{k) ■ 
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Now take the logarithm of both sides and using the ap- 
proximation log(l + x) ~ X for smaU x, we have 



N' 



log(p) « iog(D(iv)) + y: (40) 



We now approximate the sum in (j40p by an integral, so 
that 



N' 



log(P) « logiDiN)) + / (1 - 4) 



*:=1 



V + /(fc) 



Making a change of variable from k to /, gives 
log(P) « log(i?(iV)) 

(l-<5(/))-^. 



(41) 



The analysis of this equation depends upon the value of 
p and we consider separately the cases of p = 1/2 and 

5.2. The asymptotic form of the equations when p — 1/2 
We now look at the above equation when p = 1/2, and 
first determine the relation between rjk and 5^? Let 

SWp^k = Wp^k+i - Wp^k 
S log Wp^k = log Wp^k+i - log Wp^k- 



Since 



we have 



W,,k = (1 - Sk)SWp,k + Wp,k, 



logVF,,fc = log((l - 4)^^ + 1) +logW^p,fe. 
Wp^k 

Comparing with 

logW^z,fe = (1 - Vk)S\ogWp^k +logWp,fc, 
and taking Taylor expansions, for M — > oo, we have 



M 



(i-^fe)E 

m 

M 

T.(^~Sk 



(-1)'"+! fSWp^k 



w„ 



m— 1 ^ ^- 

(-1)™+! fSWp^k 



ml V W„.k 



(42) 



m—l ^ ^' 

When (SWp^k/Wp^k)^ < (SWp^k/Wp^k) it follows that 

m « Sk- 

Assuming the log poles have a normal distribution then 
6logWp,k ^ 0{1/N). For a sufficiently large network, 
when p — 0.5, we expect Sk ~ rjk for most k (the first or- 
der Taylor expansion becomes invalid near the tails of the 



normal distribution, but this contributes relatively little 
to the summation in (UHl))- The results imply that 5k is 
very close to being constant at 1/2, so that in (PT|) we have 
1- J = 1/2. 

We can then integrate the expression for P exactly. 
This allows sharp estimates of the asymptotic behaviour 
in this critical case. Integrating (|4T|) gives 



1 



log(F)«log(i^(7V)) + -lo 



/i -I- Mp^N' 



so that 



P w D{N) 



^l + Mp^N' 



In this critical case it is equally likely that we will/ will 
not have percolation paths along yi or ?/2 bonds at both 
small and large values of Accordingly, we must con- 
sider four equally likely cases of the distribution of the 
poles and zeros which could arise in any random realisation 
of the network. Thus to obtain the four possible responses 
of the network we must now consider the contribution of 
the first zero and also of the last zero. 

Case 1: First zero at the origin, last zero at iV' + 1. 
This corresponds to there being a percolation path through 
the y2 bonds. 

To determine this case we multiply P by /i to give Yi (^) 
so that 

1 

' fl + Mp^N'^ 



Mr, 



(43) 



Case 2: First zero not at the origin, last zero at N' + 
1. This corresponds to the existence of percolation paths 
through yi bonds and ?/2 bonds. 

In this case we multiply P by /i -f M^^i to give |F(^)|. 
We also use the result from the previous section that asymp- 
totically Mz^i ~ Mp.i. This then gives 



Y2{^i) « D{N)2 {^l + Mp^N')^ {p + Mp^iY' 



(44) 



Case 3: First zero at the origin, last zero at N'. Here 
there are no percolation through either set of bonds. 

To determine this case we multiply P by /x and divide 
hy n + Mz^N' to give Y . Exploiting the fact that asymp- 
totically Mp,Mi ~ Mz^N' we then have 



^3(^)^-0(^)3 



iji + Mp^N')' ifJ^ + Mp^iY' 



(45) 



Case 4: First zero not at the origin, last zero at iV'. 
This final case there exists percolation via the yi bonds 
but not through the y2 bonds. 

To determine this case we multiply P by ^ -I- M^^i and 
divide hy fi + M^.n' to give Y. Again, exploiting the fact 
that asymptotically Wp^N' ~ Wz,n' we have 



Yi{ii)^D{N), 



M, 



p,N' 



(46) 
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We know, further, from the calculations in the previous 
section that for all sufficiently large values of N 

Mp^i^l/N and Mp^jf ^ N. 

Substituting these values into the expression for Yi gives 

y,W.Z,(A,),,(ji±^)\ (47) 

The value of the constant D{N)i can be determined by 
considering the mid range of each of these expressions in 
the PLER The results of the classical Keller duality theory 
[lol | predict that each of the expressions Yi takes the same 
form in the range <C <C with 

Y,itJ-)~Vym, » = 1,2,3,4. (48) 

In the case of Yi we see that the mid-range form of the 
expression (|T7)) is given by Yi = Di^/Nji = VN^/yl/^/yi- 
This then implies that Di — yi/VN so that 

Very similar arguments lead to the following expressions 
in the other three cases: 

r2(M)«^ (iv + /i)Mi/A^ + M)s (50) 

13 (m) « ^yi r, (51) 

The four formula above give a very complete asymptotic 
description of the response of the binary network when 
p = 1/2. In particular they allow us to see the transition 
between the power-law emergent region and the percola- 
tion regions and they also describe the form of the expres- 
sions in the percolation regions. We see a clear transition 
between the emergent and the percolation regions at 

^1 = 1/iV and fi2 = N. (53) 

Hence, the number of components in the system for p = 
1/2 has a strong influence on the boundaries of the emer- 
gent region and also on the percolation response. However 
the emergent behaviour itself is independent of N. Observe 
that these frequencies correspond directly to the limiting 
pole and zero values. This gives a partial answer to the 
question, how large does N have to be to see an emergent 
response from the network. The answer is that iV has to be 
sufficiently large so that 1 /N and N are widely separated 
frequencies. 

The behaviour in the percolation regions in then given 



by the following: 







«1)^ 


a 1/2^77, 


Yi{ 


l/^l 


»1)« 




(54) 






«1). 


^/]v' 




l/^l 


»1). 


^ Va' 


(55) 






«1). 




Ysi 




»1). 


^2/iVa, 


(56) 


n( 




«1). 


2/1 


Y,{ 


l/^l 


»1). 


ayiVA. 


(57) 



We note that these percolation limits, with the strong de- 
pendence upon a/a are exactly as observed in Section [21 



5.3. The network response when p ^ 1/2 

This case differs from the case of p = 1/2 in a num- 
ber of ways and the spectral analysis is both harder and 
less complete. Firstly, rather than getting four different 
responses we expect to see only two. When p > 1/2 then 
there will (with probability one) always be conducting ca- 
pacitative percolation paths for large values of w and for 
small values of w we will not get any resistive percolation 
paths. This corresponds to the case of a first zero on the 
origin and a final zero at N' + 1. Similarly, if p < 1/2 
then we will get (with probability one) a response with 
no capacitative percolation path at high frequencies and 
resistive percolation paths at low frequencies, which corre- 
sponds to a first zero away from the origin and no final zero 
at A^'. Hence, we need only consider Case 1 and Case 4 
respectively. Secondly the values for the conductance at 
high and low frequencies are asymptotically independent 
of (sufficiently large) A^. Furthermore the formula for P 
in (|4ip involves a quadrature involving 1 — 5 which can- 
not be obtained in closed form. As a consequence we shall 
adopt a different approach for this case by combining the 
spectral calculation with that of averaging. 

6. Effective medium (averaging) calculations 

The Effective Medium Approximation (EMA) formula 
derived by an averaging method 4], gives an approxima- 
tion to the conductance of the network, and is derived by 
regarding the random distribution of the random resistors 
and capacitors as a series of perturbations of a uniform 
field identical conductors. The conductance of the effec- 
tive medium is chosen to minimise the first moment of 
the resulting perturbation matrix. It assumes an infinitely 
large number of conductances and hence corresponds to 
taking A — > oo in the previous analyses. Whilst accurate 
for p not too close to 1/2 it has limitations for p close to 
1/2, in that whilst it predicts a transition from emergent 
to percolation type behaviour, the form of this transition 
is not quite correct as p — )■ 1/2. Thus the EMA calcu- 
lations are complimentary to those derived using spectral 
methods in the previous section. In this section we will re- 
view the EMA result, and show that it is consistent with 
a PLER description of the behaviour with a power law 
which we explicitly derive. We we then make a (somewhat 
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speculative, but consistent) extension of the EM A calcu- 
lation to include the effects of finite network size N. We 
see presently that if TV > A^* = \p-l/2\-'^ then the EMA 
formula gives a good approximation to the resulting con- 
ductance and we will obtain a more general formula which 
is effective for all |p — 1/2| and 1/VN sufficiently small. 

6.1. Infinite networks 
6.1.1. Overview 

If the conductances are yi and ?/2 are in respective pro- 
ortion 1 — p and p then, from the 'classical' EMA result, 
, the effective medium conductance Y for a very large 
{N — )• oo) square two-dimensional lattice is given by the 
solution of the quadratic equation 



y-yi 



p 



Y-y2 



0. 



Y + yiJ '\Y + y2^ 
Rearranging this formula we have 

Y^ + {l-2p){y2-yi)Y-y,y2^0, 

so that if 

e^{l-2p), 9 = Y/^yiy2, n^y2/yi, 
we have 



(58) 



Now set 



7 = log (6*), and = log(/Lt) 
it follows immediately that 

sinh(7) = — e sinh(zy/2). 



(60) 



6.1.2. Emergent power laws 

Suppose firstly that /i is reaZ, so that we are modelling 
a R-R network, and that /i is close to unity, so that v, and 
hence 7 are both not large. Then we may linearise (|60p 
and to leading order we have 7 — —ev 12. Thus in this case 
log(F/y^yIyJ) = elog(?/i/?/2)/2, and rearranging this gives 
the elegant power law identity 



(61) 



This is fully consistent with the duality result that 

y {V\,V2)Y {y2,V2) = V'x^'^^ vlvlvh"^^ = V\V2- 

In the C-R network case, /i = iu)CR is pure imaginary. 
We set ^ = ir] where r] — ujCR is now assumed to be close 
to unity and take /3 — \og{r]) to be close to zero. It then 
follows that log(/x) — i-K/2 + P so that 

sinh(7) -esinh(i7r/4-|-/3/2) = -e{i + (3/2 + 0{l3'^))/V2. 

If /3 = then 7 — iOo where sin(0o) — — e/v^- Linearising 
about this solution we then have, to leading order. 



7 = i«o - 



2%/2cos(6'o) 



/3 + 0(/32) = tOo -f Alog(??) + 0(/32). 



Thus, to leading order 



|r| = y/ujC/R I exp(7)| = y/ijC/R {loCR)^ = Kuj". 

This gives precisely the power law observed in the results 
presented in Section [5] in a frequency range centred around 
the region for which uCR = 0(1) and 



a(p) = - + A = 

2 2 2V2 



l-eV2 



(62) 



Note that ii p = 1/2, so that e = 0, then we see power 
law behaviour with exponent 1/2 in all cases, over the 
entire range. This contrasts with the results of the last 
section, in which we see percolation type behaviour for 
(say) u! > N, but this also emphasises the fact that the 
EMA calculation only applies for A^ = cc in this case. 
Note further that a(0) = 0,a(l/2) = 1/2, a(l) = 1 and 
that a{p) « p for all < p < 1. 

6.1.3. Percolation behaviour 

It is well known that the EMA approximation for p 7^ 
1/2 exhibits percolation behaviour which we summarise 
here. If we consider the case of large positive v then the 
asymptotic form of the solution depends upon the sign of 
e. If e > 0, which corresponds to p < 1/2 then for large a 
the equation ([SO)) reduces to e^'^ — ee"/"^ so that we have 
the percolation behaviour given by: 



l/e - eVTZ, r = yi/e 



(63) 



In contrast, ife<0, (p>l/2) then for large v the equa- 
tion ((60|) simplifies to: e'' — {—e)e^/'^ so that wc have the 
percolation behaviour given by: 



(64) 



Similar results for v large and negative follow from duality 
arguments. The (frequency) limits of the emergent region 
can be estimated by finding when the power law behaviour 
of (l6lT) overlaps with the percolation type behaviour. This 
leads to the following estimates for the values fii < ^ < 112 
over which we expect to see power-law emergent behaviour 



e > 
e < 



Ml 
Ml 



i/(i-p) 



(65) 



Important Note These results predict that as e — >■ 
the percolation amplitude scale as |e|^^. This is not 
quite what is observed in practice. In contrast, empirical 
calculations described in (for example) @ , imply instead 
a scaling law of the from jel^^ "^. Thus the EMA is not fully 
accurate in this limit. 

6.2. Large, but finite, networks 

We now give a more speculative calculation which at- 
tempts to combine the EMA estimate with finite size ef- 
fects and the spectral calculations of the previous section. 



16 



for Case 1 and Case 4. Our starting point is the spectrally 
derived formula Yi which has percolation limits pro- 
portional to ?/2 when fj, is large. Casting Yi in terms of yi 
and 2/2 we have 



N {(1 + 1/N) 
(1 + fi/N) 



(1 + l/N^i) 

(1 + t^/N) 

(1 + l/N^l) ■ 



It follows that 



(1 + l/N^l)Y^ - yiy2(l + m/A^) = 0. 

If, as above, we set 6 = Y/ -^2/12/2 we have, after some 
manipulation, that this formula has the symmetric form 



11/^ 9 

0^ N [j^ Jl 



(66) 



Similarly, the spectrally derived formula Y4 in (j52p , which 
has percolation limits proportional to yi for /i large, takes 
the symmetric form 



1 



1 



1 



fi0 



(67) 



These expressions are both very similar in form to the 
result ([5^ of the EMA calculation. We conjecture that 
a more general expression can be obtained by combining 
them into the following two fomulae which agree with each 
in the limits of e = and N — 00 and which include both 
the effects of component proportion and network size and 
which respectively have percolation limits proportional to 
yi and 7/2: 



1 



1 / 1 



1 



(68) 



and 



(69) 



We observe that each of and (IM)) is self-dual un- 
der the map /i — > 9 — > 1/9. Similarly the symme- 
try ^ — > 1/^, e — > — e maps (^5]) to ([55)1 and vice versa. 
We now proceed to show that (j68l69l) have the correct 
asymptotic form of solution and give numerical evidence 
for their validity in Section [71 We firstly consider the per- 
colation limits of ([55[) and ([55[) . Motivated by the analy- 
sis in the previous subsection we consider solutions of the 
form 9 = PyTp, so that Y = f3y2, and 9 = f3/y/JI, so that 
Y = Pyi, in the two cases of fi large and fi small. If fi 
is large and e > then ([55)1 has a solution with percola- 
tion limit proportional to, and in phase with, j/i, so that 
= PI if /3 satisfies the quadratic equation 



^+e/3-l = 0. 



(70) 



If /z is small then we have the reciprocal solution given by 
the map /3 ^ 1//3. Note that if 6* = /?/^then |r| = /3|?/i| 
and hence the dynamic range is given by 



fl2 

Y = 0^ where ^ + e/? 
It is immediate that /3 is given by 



N 



(71) 



(72) 



where the positive sign for the square root term is taken 
to ensure that /3 > so that the response is in phase with 
yi. This expression takes two different forms depending 
on whether (i) N or (ii) N :$> e~^. In the first 

case the network behaves in a similar way to one with 
p = 1/2 and we have (3 ^ VTV- In the second case we have 
behaviour similar to = 00 with /3 ^ 1/e. This is in exact 
correspondence with the calculations of the dynamic range 
reported in Section [2| Similarly, if fi is large and e < 
then the equation (j69p has a solution with percolation limit 
proportional to ,and in phase with, y2, so that 9 — 
if P satisfies the quadratic equation 



P' 



0. 



(73) 



7. Comparison of the asymptotic and numerical 
results 

We now give two sets of calculations for finite networks. 
The first tests the validity of the spectral calculation in 
Section[51for the case oip =1/2. The second the validity of 
the amalgamated spectral and averaging based calculation 
in Section [6| 

7.1. Spectral based calculations for p = 1/2 

We firstly consider the case oip = 1/2 (e = 0) for the C- 
R network with complex /i = iujCR. We compare the ab- 
solute values of the four asymptotic formute ([49I50I51I52P 
obtained by using the spectral method with the numerical 
calculations of the absolute network conductance |F| with 
C = InF and R = IkQ, as a function of lo for four dif- 
ferent configurations of the system, with different percola- 
tion paths for low and high frequencies. The results of this 
comparison are shown in Figure [ijl in which we plot the 
numerical calculations together with the asymptotic for- 
mulae for a range of values of N given by = S'(S' — 1) with 
S = 10, 20, 50, 100. We can see from this that the predic- 
tions of the asymptotic formulae (I49I50I51I52P fit perfectly 
with the results of the numerical computations over all of 
the values of A^ considered. Indeed they agree both in the 
(square-root) power law emergent region and in the four 
possible percolation regions. The results and the asymp- 
totic formula clearly demonstrate the effect of the network 
size in these cases. 

As a separate calculation, we look at the results ob- 
tained for the R-R network with conductances 1/R and 
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/x/i? with real /i and with R as above. Again we compare 
the predictions of the asymptotic formulas (|49I50I51I52I) 
with the numerical computations of |y | in this case. Again 
we see an excellent agreement in all cases, as shown in Fig- 
ure [131 

7.2. Combined averaging and spectral based calculations 
for general p 

We now consider the responses described by the pair of 
quadratic equations (j68l69p for Y obtained by combining 
the averaging and spectral calculations. We compare these 
with numerical results for the C-R networks described in 
the previous sub-section. In each case of p we take the 
equation for which the corresponding solution in the per- 
colation regime is physically correct. 

As a first computation we take p ~ 0.4, so that e = 
0.2 > 0, and consider a C-R network with the same val- 
ues of C,R and taking S so that N = 90,9900. For an 
infinitely large network we expect to see resitive type per- 
colation behaviour (with probability one) for both small 
and large frequencies. We compute from (|55)) and 

compare these values with the results of computations of 
|y(cLi)| from a large realisations of the network in Figure 
[TBI The results from this computation are interesting. We 
see in the case of = 90, that there is quite a large sta- 
tistical range in the calculations in this case. In the case 
of A^ = 9900 the statistical range is much smaller. The 
results of the predictions of |y| from ([68l) closely follow 
the mid-range of the computations for A^ = 90. In the 
case of A^ = 9900 there continues to be good agreement, 
however, as expected from the EMA results, the maximum 
and minimum values of |K| are slightly underestimated by 
([68)) . In both cases the PLER is very well approximated 
by the computations from (j68p . The results for computa- 
tions of the R-R network are very similar and we do not 
include them here. 

As a second computation we take p = 0.6, so that 
e = —0.2 < 0. For an infinitely large C-R network we 
expect to see reactive type percolation behaviour (with 
probability one) for both small and large frequencies. We 
present the results of computing |F(a;)| from the quadratic 
equation (j69p compared with computations from a number 
of realisations of the network when A^ = 90, 9900, in Figure 
1171 In this computation we again see a greater statistical 
range when A^ = 90 than when A^ — 9900. Indeed in the 
case of A^ = 90 a small number of the realisations show 
resistive percolation behaviour rather than reactive. This 
is not seen in the calculations for A^ — 9900. In both cases 
the results of the calculations from ([69)) closely match the 
computations over the whole range. 



Figure 14: (color online) Comparison of the asymptotic formulee ob- 
tained from the spectral derivation with the numerical computations 
for the C — R network over many runs, with p = 1/2 and network 
sizes sizes S = 10, 20, 50, 100, N = S{S - 1). 



8. Discussion 

We have attempted to answer questions surrounding 
emergent behaviour, determining what causes it, finding 
the range of parameters over which it applies and address- 
ing the question of which aspects of a complex system 
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Figure 15: (color online) Comparison of the asymptotic formulas 
with the numerical computations for the R — R network over many 
runs, with p = 1/2 and network sizes sizes S = 10, 20, 50, 100, N = 
S{S~1). 
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Figure 16: (color online) Comparison of the numerically obtained 
solutions of the quadratic equation II68I I with the numerical compu- 
tations for many realisations of the the C — R network, with p = 0.4 
and network sizes sizes S = 10, 100, N = S{S — 1). 



influence the emergent behaviour. Using large binary net- 
works we have shown how power law emergence can be di- 
rectly related to the statistical regularity of the spectrum 
of the matrices associated with the network and hence can 
be studied by combing spectral and averaging methods. 
In particular we have studied the effects of network size, 
and the variation from criticality on the observed power 
law behaviour of these systems. We have shown how the 
response of the networks depends strongly upon p and less 
strongly on the network size N, except at p = 1/2 exactly, 
where the dynamic range has been found to scale in direct 
proportion to N. When p = 1/2 we analysed how the net- 
work response is described in terms of poles and zeros of 
the conductance and can be determined from distribution 
of these values, making use of numerically observed statis- 
tical patterns of these. This has revealed four asymptotic 
formulae, corresponding to the four qualitatively different 
emergent responses that can arise when p — 1/2 and these 
show very precisely the effects of the (finite) network size 
N. The case of p = 1/2 is very complete asymptotically 
and shows particularly good agreement with the numeri- 
cal computations, which is remarkable given the number 
of approximations made. An important open question is 
to now rigorously establish the observed statistical results 
of the spectrum in this case, for example to show rigor- 
ously that fip^N' ~ N. When p ^ 1/2 the analysis is less 
complete. It is interesting, however, that the results of the 
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Figure 17: (color online) Comparison of the numerically obtained 
solutions of the quadratic equation II69I I with the numerical compu- 
tations of 1 5^(1^)1 for many realisations of the C — R network with 
p = 0.6 and network sizes S = 10, 100, = S(S - 1). 

averaging based EMA calculation can be combined with 
those of the spectral computation in a consistent manner 
to the case of finite iV, leading to predictions (|68I69L of 
the conductance and its dynamic range which is in good 
qualitative agreement with what is observed. However a 
limitation of this analysis remains the lack of precision of 
the estimation of the power law scaling of the magnitude 
of the percolation response as p ^ 1/2. We conclude that 
combining both the spectral based and the averaging based 
methods lead to useful asymptotic formulas with excellent 
numerical support, and establishing these more rigorously 
is an interesting area of further study. 
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